#include<iostream>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX、AVX2
#include<windows.h>
#include<stdlib.h>
using namespace std;
int N = 2000;

//生成不会出错的测试用例
void m_reset(float** m, int n)
{
    for (int i = 0; i < n; i++)
    {
        // 下三角全部设置为0
        for (int j = 0; j < i; j++)
            m[i][j] = 0;
        //对角线设置为1.0
        m[i][i] = 1.0;
        //上三角设置为随机数
        for (int j = i + 1; j < n; j++)
            m[i][j] = rand() % 100;
    }

    //每一行加上比自己下标小的行
    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
            for (int j = 0; j < n; j++)
                m[i][j] += m[k][j];
}

//并行
void parallel_solve(float** m, int n)
{
    for (int k = 0; k < n; k++)
    {
        //1、先是第k行除以m[k][k]
        //Broadcast single-precision (32-bit)
        //floating-point value a to all elements of dst.
        __m256 vt = _mm256_set1_ps(m[k][k]);
        int j = k + 1;

        //开头对齐
        for (; j < n; j++) {
            if (j % 8 == 0) break;
            m[k][j] = m[k][j] / m[k][k];
        }

        for (; j + 8 <= n; j = j + 8)
        {
            //对齐模式
            __m256 va = _mm256_load_ps(&m[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_store_ps(&m[k][j], va);
        }
        //处理剩余的
        if (j < n)
        {
            for (; j < n; j++)
                m[k][j] = m[k][j] / m[k][k];
        }


        //1、处理m[k][k]
        m[k][k] = 1.0;

        //消除第k列行下标大于k的元素
        for (int i = k + 1; i < n; i++)
        {
            __m256 vaik = _mm256_set1_ps(m[i][k]);
            int j = k + 1;
            //处理开头部分
            for (; j < n; j++)
            {
                if (j % 8 == 0) break;
                m[i][j] = m[i][j] - m[i][k] * m[k][j];
            }

            //向量计算
            for (; j + 8 <= n; j = j + 8)
            {
                __m256 vakj = _mm256_load_ps(&m[k][j]);
                __m256 vaij = _mm256_load_ps(&m[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_store_ps(&m[i][j], vaij);
            }
            //处理剩余的
            if (j < n)
            {
                for (; j < n; j++)
                    m[i][j] = m[i][j] - m[k][j] * m[i][k];
            }
            m[i][k] = 0;
        }

    }
}

int main()
{
    float** matrix_parallel;
    matrix_parallel = new float* [N];
    for (int i = 0; i < N; i++)
    {
        matrix_parallel[i] = new float[N];
    }

    int size = 32;//数据集增长速度
    long long head1, tail1, freq1;//记录串行时间
    for (int n = size; n < N; n += size)
    {
        m_reset(matrix_parallel, n);

        //并行
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq1);
        QueryPerformanceCounter((LARGE_INTEGER*)&head1);
        QueryPerformanceCounter((LARGE_INTEGER*)&tail1);
        int count1 = 0;
        while ((tail1 - head1) < 0.1 * freq1) {
            parallel_solve(matrix_parallel, n);
            QueryPerformanceCounter((LARGE_INTEGER*)&tail1);
            count1++;
        }
        float time1 = (tail1 - head1) * 1000.0 / (freq1 * count1);

        //并行

        cout << "问题规模:" << n << endl;
        cout << "对齐并行平均时间" << time1 << "ms";
        cout << endl;

        if (n > 320) size = 100;

    }

}

